

%Run programs
clc
clear all
close all
% 
%Initial steady-state
run Parameters
run BGP_growth

ss_c=C;
ss_y=Y;
ss_yT=YT;
ss_yNT=YNT;
ss_ha=Ha;
ss_hr=H_r;
ss_g=g(length(g));
ss_TT=Tnew;
ss_TNT=TNT;
ss_pi=Pi;
ss_piT=PiT;
ss_piNT=PiNT;
ss_p=P;
ss_pT=PT;
ss_pNT=PNT;
ss_PishareT=pi_shareT;
ss_PishareNT=pi_shareNT;
ss_w=w;
ss_vinnovator=vinnovator;
ss_jinnovator=Jinnovator;
ss_v=vadopter;
ss_j=Jadopter;
ss_vinnov=vinnov;
ss_chi = chi;
ss_YL=Y_L;
ss_eps=eps;
ss_Z=ZZ;
ss_AT=A_Z;
ss_rp=roy;
ss_IM=LHS;
ss_EX=RHS;
ss_IMT=IM_T;
ss_EXT=EX_T;
ss_IMNT=IM_NT;
ss_EXNT=EX_NT;
ss_Yw=Yw;
ss_Deltaeps=Delta_eps;
ss_rhochi = rho_chi; 

save Initial_values.mat ss_rhochi ss_chi ss_Yw Q ss_TNT ss_y ss_yT ss_yNT ss_pT ss_pNT ss_IMT ss_EXT ss_IMNT ss_EXNT ss_Deltaeps ss_c ss_y ss_ha ss_hr ss_g ss_TT ss_pi ss_piT ss_piNT ss_p ss_PishareT ss_PishareNT ss_w ss_vinnovator ss_jinnovator ss_v ss_j ss_vinnov ss_YL ss_eps ss_Z ss_AT ss_rp ss_EX ss_IM;
save Q.mat Q
save epsbar.mat epsbar
save lam.mat lam


%%Counterfactual steady-state
run Parameters_counterf
run BGP_growth_counterf


ss1_c=C;
ss1_y=Y;
ss1_yT=YT;
ss1_yNT=YNT;
ss1_ha=Ha;
ss1_hr=H_r;
ss1_g=g(length(g));
ss1_TT=Tnew;
ss1_TNT=TNT;
ss1_pi=Pi;
ss1_piT=PiT;
ss1_piNT=PiNT;
ss1_p=P;
ss1_pT=PT;
ss1_pNT=PNT;
ss1_PishareT=pi_shareT;
ss1_PishareNT=pi_shareNT;
ss1_w=w;
ss1_vinnovator=vinnovator;
ss1_chi = chi;
ss1_jinnovator=Jinnovator;
ss1_v=vadopter;
ss1_j=Jadopter;
ss1_vinnov=vinnov;
ss1_YL=Y_L;
ss1_eps=eps;
ss1_Z=ZZ;
ss1_AT=A_Z;
ss1_rp=roy;
ss1_IM=LHS;
ss1_EX=RHS;
ss1_IMT=IM_T;
ss1_EXT=EX_T;
ss1_IMNT=IM_NT;
ss1_EXNT=EX_NT;
ss1_Yw=Yw;
ss1_Deltaeps=Delta_eps;
ss1_rhochi = rho_chi; 

save Final_values.mat ss1_rhochi ss1_chi ss1_Yw ss1_TNT ss1_y ss1_yT ss1_yNT ss1_pT ss1_pNT ss1_IMT ss1_EXT ss1_IMNT ss1_EXNT ss1_Deltaeps ss1_c ss1_y ss1_ha ss1_hr ss1_g ss1_TT ss1_pi ss1_piT ss1_piNT ss1_p ss1_PishareT ss1_PishareNT ss1_w ss1_vinnovator ss1_jinnovator ss1_v ss1_j ss1_vinnov ss1_YL ss1_eps ss1_Z ss1_AT ss1_rp ss1_EX ss1_IM;
save epsbar.mat epsbar

 clear all

  dynare 'Dynamic_GFT_exponents.mod' savemacro 
  save DynamicGFT.mat

%%Compute statistics

load Initial_values.mat
load Final_values.mat

dYL=log(ss1_YL./ss_YL)'.*100+(ss1_g-ss_g)*100/4
dC=log((ss1_c./ss1_p')./(ss_c./ss_p')).*100+(ss1_g-ss_g)*100/4./(1-bet)
g0=ss_g
g1=ss1_g
dhts=log(diag(ss1_PishareT)./diag(ss_PishareT)).*100
dimport_share=log((1-diag(ss1_PishareT))./(1-diag(ss_PishareT))).*100
dhr=log(ss1_hr./ss_hr)'.*100;
dha_total=log(sum(ss1_ha')./sum(ss_ha'))'.*100;
dhr_y=log((ss1_hr'./ss1_y)./(ss_hr'./ss_y)).*100;
dT=log((ss1_TT./ss1_TT(3))./(ss_TT./ss_TT(3))).*100;
dha_y=log((sum(ss1_ha')./ss1_y')./(sum(ss_ha')./ss_y'))'.*100;
Tss0=ss_TT./ss_TT(M);
Tss1=ss1_TT./ss1_TT(M);
%Dynamic gains from trade

 load DynamicGFT.mat
% 
TT=150;
betg=bet; 
kk=[1:TT];

for i=1:numit

% Uold1(i)=sum((betg.^([1:TT]-1))'.*log(exp(resultsC1(i,1))))+sum((betg.^([1:TT]-1))'.*(cumsum(ones(TT,1).*(1+resultsgc1(i,1)))));
% Unew1(i)=sum((betg.^([1:TT]-1))'.*log(exp(resultsC1(i,1:TT)')))+sum((betg.^([1:TT]-1))'.*(cumsum(ones(TT,1).*((1+resultsgc1(i,1:TT)'./(1))))));
% 



Uold1(i)=sum((betg.^([1:TT]).*log((1+resultsgc1(i,1)).^[1:TT])));
prodgc1 = cumprod((1+resultsgc1(i,1:TT)));
Unew1(i)=sum((betg.^([1:TT]).*log(prodgc1)));



%GFT1(i)=[(exp((1-betg)*(Unew1(i)-Uold1(i))))-1]*100;


% Uold3(i)=sum((betg.^([1:TT]-1))'.*log(exp(resultsC3(i,1))))+sum((betg.^([1:TT]-1))'.*(cumsum(ones(TT,1).*(1+resultsgc3(i,1)))));
% Unew3(i)=sum((betg.^([1:TT]-1))'.*log(exp(resultsC3(i,1:TT)')))+sum((betg.^([1:TT]-1))'.*(cumsum(ones(TT,1).*((1+resultsgc3(i,1:TT)'./(1))))));


Uold3(i)=sum((betg.^([1:TT]).*log((1+resultsgc3(i,1)).^[1:TT])));
prodgc3 = cumprod((1+resultsgc3(i,1:TT)));
Unew3(i)=sum((betg.^([1:TT]).*log(prodgc3)));


Uold2(i)=sum((betg.^([1:TT]).*log((1+resultsgc2(i,1)).^[1:TT])));
prodgc2 = cumprod((1+resultsgc2(i,1:TT)));
Unew2(i)=sum((betg.^([1:TT]).*log(prodgc2)));



GFT3(i)=[exp((1-betg)*((log(exp(resultsC3(i,TT))*(1+resultsgc3(i,TT)))-log(exp(resultsC3(i,1))*(1+resultsgc3(i,1))))/(1-betg)+Unew3(i)-Uold3(i)))-1]*100;
GFT1(i)=[exp((1-betg)*((log(exp(resultsC1(i,TT))*(1+resultsgc1(i,TT)))-log(exp(resultsC1(i,1))*(1+resultsgc1(i,1))))/(1-betg)+Unew1(i)-Uold1(i)))-1]*100;
GFT2(i)=[exp((1-betg)*((log(exp(resultsC2(i,TT))*(1+resultsgc2(i,TT)))-log(exp(resultsC2(i,1))*(1+resultsgc2(i,1))))/(1-betg)+Unew2(i)-Uold2(i)))-1]*100;


C1counterf(i,:)=log(exp(resultsC1(i,TT))*(1+resultsgc1(i,TT))^(sig-1))-log(exp(resultsC1(i,1))*(1+resultsgc1(i,1))^(sig-1))+log((prodgc1))-log((1+resultsgc1(i,1)).^[1:TT]);
C3counterf(i,:)=log(exp(resultsC3(i,TT))*(1+resultsgc3(i,TT))^(sig-1))-log(exp(resultsC3(i,1))*(1+resultsgc3(i,1))^(sig-1))+log((prodgc3))-log((1+resultsgc3(i,1)).^[1:TT]);
C2counterf(i,:)=log(exp(resultsC2(i,TT))*(1+resultsgc2(i,TT))^(sig-1))-log(exp(resultsC2(i,1))*(1+resultsgc2(i,1))^(sig-1))+log((prodgc2))-log((1+resultsgc2(i,1)).^[1:TT]);

end


save Results_NN.mat


load Results_NN.
p=1021;

close all
hold on
plot([0 0 0 0 0 0 0 0  0 C1counterf(p,1:150)])
plot([0 0 0 0 0 0 0 0  0 C3counterf(p,1:150)])
%plot([0 0 0 0 0 0 0 0  0 C2counterf(340,1:50)])
hold off


% 
% for i=1:11
% for j=1:length(GFT1)
% GFT(i,j) = rho(i)*GFT1(j)+(1-rho(i))*GFT3(j);
% GFT1_n(i,j) = GFT1(j);
% GFT3_n(i,j) = GFT3(j);
% end
% end

